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Nonuniform neutron-rich matter present in both core-collapse supernovae and neutron-star crusts 
is described in terms of a semiclassical model that reproduces nuclear-matter properties and includes 
long-range Coulomb interactions. The neutron-neutron correlation function and the correspond- 
ing static structure factor are calculated from molecular dynamics simulations involving 40,000 to 
100,000 nucleons. The static structure factor describes coherent neutrino scattering which is ex- 
pected to dominate the neutrino opacity. At low momentum transfers the static structure factor is 
found to be small because of ion screening. In contrast, at intermediate momentum transfers the 
static structure factor displays a large peak due to coherent scattering from all the neutrons in a 
cluster. This peak moves to higher momentum transfers and decreases in amplitude as the density 
increases. A large static structure factor at zero momentum transfer, indicative of large density 
fluctuations during a first-order phase transition, may increase the neutrino opacity. However, no 
evidence of such an increase has been found. Therefore, it is unlikely that the system undergoes a 
simple first-order phase transition. Further, to compare our results to more conventional approaches, 
a cluster algorithm is introduced to determine the composition of the clusters in our simulations. 
Neutrino opacities are then calculated within a single heavy nucleus approximation as is done in most 
current supernova simulations. It is found that corrections to the single heavy nucleus approxima- 
tion first appear at a density of the order of 10 13 g/cm 3 and increase rapidly with increasing density. 
Thus, neutrino opacities are overestimated in the single heavy nucleus approximation relative to the 
complete molecular dynamics simulations. 

PACS numbers: 



I. INTRODUCTION 



The description of nuclear matter at subnuclear den- 
sities is an important and general problem. Attrac- 
tive short-range strong interactions correlate nucleons 
into nuclei. However, nuclear sizes are limited by long- 
range repulsive Coulomb interactions and thermal exci- 
tations. This competition between attraction and repul- 
sion produces multifragmentation in heavy ion collisions; 
the breaking of the system into intermediate sized frag- 
ments [l|, |2|, |3| . I n astrophysics this competition produces 
a variety of complex phenomena. At densities consider- 
ably lower than normal nuclear matter saturation den- 
sity, the system may be described as a collection of nearly 
free nucleons and nuclei in nuclear statistical equilibrium 
(NSE), while at normal nuclear matter saturation den- 
sity and above the system is expected to become uniform. 
In between these regimes pasta phases may develop with 
nucleons clustered into subtle and complex shapes 0, 0| . 
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This pasta phase may be present in the inner crust of 
neutron stars and in core collapse supernovae. Unfortu- 
nately, low-density NSE models and high-density models 
of nuclear matter are often incompatible. As the den- 
sity increases, it is difficult to account for the strong in- 
teractions between nuclei in NSE models. Likewise, as 
the density decreases and the uniform system becomes 
unstable against fragmentation, uniform models fail to 
describe cluster formation. 

We wish to study the different phases and properties 
of the system as a function of density. This is essential 
in the simulations of core-collapse supernovae as they in- 
volve a tremendous range of densities and temperatures. 
Several semiclassical simulation techniques, often devel- 
oped for heavy ion collisions, may be used to describe the 
system over this large density range. Indeed, Watanabe 
and collaborators have used quantum molecular dynam- 
ics to describe the structure of the pasta @|. These sim- 
ulations should reduce to isolated nuclei at low densities 
and to uniform matter at high densities. In principle a 
first order phase transition could have a two phase coexis- 
tence region. Large density fluctuations at the transition 
could greatly increase the neutrino opacity 0. In the 
present paper, we search for regions with large density 
fluctuations using molecular dynamics simulations. 

The main focus of the present paper is coherent neu- 
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trino scattering, an essential tool for computing neutrino 
mean-free paths in supernovae and to determine how neu- 
trinos are initially trapped. Furthermore, coherent neu- 
trino scattering — by representing long-range "classical" 
physics — may provide insights into how the clustering 
evolves with density in a model-independent way. In a 
previous paper a simple Monte-Carlo simulation model 
involving 4,000 particles was developed and first results 
for neutrino mean- free paths were presented ||. In this 
paper results are presented for larger simulations involv- 
ing up to 100,000 nucleons using molecular dynamics. 
This large number of nucleons is required to accommo- 
date long wavelength neutrinos. For example, the wave- 
length of a 10 MeV neutrino is approximately 120 fm. At 
a baryon density of 0.05 fm -3 , a simulation volume of one 
neutrino wavelength per side contains close to 100,000 
nucleons. 

Having generated a variety of observables, our micro- 
scopic results are then compared to those generated from 
a macroscopic cluster model. Such macroscopic mod- 
els describe the system as a collection of free nucleons 
plus a single species of a heavy nucleus and are presently 
used in most supernovae simulations. By comparing the 
two approaches, we gain insight into the strengths and 
limitations of the macroscopic cluster models. There is 
a duality between microscopic descriptions of the sys- 
tem in terms of nucleon coordinates and macroscopic de- 
scriptions in terms of effective nuclear degrees of free- 
dom. Thus, it is interesting to learn when does a neu- 
trino scatter coherently from a nucleus and when does 
it scatter from an individual nucleon? At the Jefferson 
Laboratory a similar question is being posed: when does 
a photon couple coherently to a full hadron and when 
to an individual quark? The quark/hadron duality has 
provided insight on how descriptions in terms of hadron 
degrees of freedom can be equivalent to descriptions in 
terms of quark coordinates |9( . Here we are interested in 
nucleon/nuclear duality: how can nuclear models incor- 
porate the main features of microscopic nucleon descrip- 
tions? 

The manuscript has been organized as follows. In 
Sec.[n]the simple semiclassical formalism is introduced as 
well as details of the molecular dynamics simulations. In 
Sec. 11111 we review the formalism for neutrino scattering 
and relate it to the static structure factor. Simulation 
results are presented in Sec. lIVI including the calculation 
of neutrino mean-free paths using nucleon coordinates. 
Section [3 presents a simple cluster model that compares 
these results to more conventional approaches using nu- 
clear coordinates. Finally, conclusions and future direc- 
tions are presented in Sec. IVII 



II. FORMALISM 

In this section we review our semiclassical model that 
while simple, contains the essential physics of competing 
interactions consisting of a short-range nuclear attraction 



and a long-range Coulomb repulsion. The impossibility 
to simultaneously minimize all elementary interactions is 
known in condensed-matter circles as frustration. The 
complex physics of frustration, along with many other 
details of the model, may be found in Ref . |g ■ Here only 
a brief review of the most essential features of the model 
is presented. We model a charge-neutral system of elec- 
trons, protons, and neutrons. The electrons are assumed 
to be noninteracting and thus are described as a degen- 
erate free Fermi gas at a number density identical to that 
of the protons (i.e., p e — p p ). The nucleons, on the other 
hand, interact classically via a nuclear-plus-Coulomb po- 
tential. However, the use of an effective temperature 
and effective interactions are used to simulate effects as- 
sociated with quantum zero-point motion. More elabo- 
rate models are currently under construction and these 
will be presented in future contributions. While simple, 
the model displays the essential physics of frustration, 
namely, nucleons clustering into pasta but the size of the 
clusters limited by the Coulomb repulsion, in a trans- 
parent form. Moreover, one may study the evolution of 
the system through the low density, pasta, and high den- 
sity phases within a single microscopic model. Finally, 
the model facilitates simulations with a large numbers 
of particles, a feature that is essential to estimate and 
control finite-size effects and, as alluded earlier, to reli- 
ably study the response of the system to long wavelength 
neutrinos. 

The total potential Vtot energy of the system consists 
of a sum of two-body interactions 



Vtot = 



i<j 



V(i,j) 



(1) 



where the "elementary" two-body interaction is given as 
follows: 



V(i,j) = ae 



b + cT z (i)r z (j) 



,/2A 



V c (i,j) 



(2) 

Here the distance between the particles is denoted by 
Tij = |r,- — Yj\ and r z represents the nucleon isospin 
projection (r 2 = +1 for protons and t 2 = —1 for neu- 
trons). The two-body interaction contains the charac- 
teristic intermediate-range attraction and short-range re- 
pulsion of the nucleon-nucleon force. Further, an isospin 
dependence has been incorporated in the potential to en- 
sure that while pure neutron matter is unbound, sym- 
metric nuclear matter is appropriately bound. Indeed, 
the four model parameters (a, b, c, and A) introduced in 
Eq. J2J have been adjusted in Ref. @ to reproduce the 
following bulk properties: a) the saturation density and 
binding energy per nucleon of symmetric nuclear matter, 
b) (a reasonable value for) the binding energy per nucleon 
of neutron matter at saturation density, and c) (approx- 
imate values for the) binding energy of a few selected 
finite nuclei. All these properties were computed via a 
classical Monte Carlo simulation with the temperature 
arbitrarily fixed at 1 MeV. The parameter set employed 
in all previous and present calculations is displayed in 
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Table [I] Finally — and critical for pasta formation — a 
screened Coulomb interaction of the following form is in- 
cluded: 



V c (i,j) = — e- r ^ X T p (i)r p (j) , 



(3) 



where t p = (l+r z )/2 and A is the screening length that 
results from the slight polarization of the electron gas. 
The relativistic Thomas-Fermi screening length is given 
by 



A = — I k F \/kl + ml 



-1/2 



(4) 



where m e is the electron mass and the electron Fermi 
momentum has been defined by k-p = (37r 2 p e ) 1 / 3 [IcLllH. 
Unfortunately, while the screening length A defined above 
is smaller than the length L of our simulation box, it is 
not significantly smaller. Hence, following a prescription 
introduced in Ref. §| in an effort to control finite-size 
effects, the value of the screening length is arbitrarily 
decreased to A = 10 fm. 

The simulations are carried out with both a fixed num- 
ber of particles A and a fixed density p. The simulation 
volume is then simply given by V — A/ p. To minimize 
finite-size effects periodic boundary conditions are used, 
so that the distance is calculated from the x, y, and 
z coordinates of the ith and j'th particles as follows: 



Xi 



+ [Vi - Vj} 2 + [«i ~ z ] 



(5) 



where the periodic distance, for a cubic box of side L 
V 1 / 3 , is given by 



[I] = Min(|Z|,L- 



(0) 



In our earlier work y| properties of the pasta were ob- 
tained from a partition function that was calculated using 
Monte Carlo integration. In the present paper molec- 
ular dynamics is used to simulate the system. There 
are a few advantages in using molecular dynamics over 
a partition function. First, larger systems are allowed 
to be simulated due to advances in both software and 
in hardware (see appendix). Second, one is not limited 
to compute the static structure factor of the system as 
now the full dynamic response is available. One expects 
the neutron-rich pasta to display interesting low-energy 
collective excitations — such as Pygmy resonances — that 
may be efficiently excited by low-energy neutrinos. These 
low-energy modes of the pasta are currently under inves- 
tigation. 

To carry out molecular dynamics simulations the tra- 
jectories of all of the particles in the system are de- 
termined by simply integrating Newton's laws of mo- 
tion, albeit for a very large number of particles (up to 
100,000 in the present case) using the velocity- Verlet 
algorithm [l2L fl3l IbH ]. To start the simulations, ini- 
tial positions and velocities must be specified for all the 



TABLE I: Model parameters used in the calculations. 



a 


b 


c 


A 


110 MeV 


-26 MeV 


24 MeV 


1.25 fm 2 



particles in the system. The initial positions are ran- 
domly and uniformly distributed throughout the simu- 
lation volume while the initial velocities are distributed 
according to a Boltzmann distribution at temperature T. 
As the velocity- Verlet is an energy — not temperature — 
conserving algorithm, kinetic and potential energy con- 
tinuously transformed into each other. To prevent these 
temperature fluctuations, the velocities of all the parti- 
cles are periodically rescaled to ensure that the average 
kinetic energy per particle remains fixed (3/2) /«bT. 

In summary, a classical system has been constructed 
with a total potential energy given as a sum of two- 
body, momentum-independent interactions as indicated 
in Eq. (J2J). Expectation values of any observable of in- 
terest may be calculated as a suitable time average using 
particle trajectories generated from molecular dynamics 
simulations. 



III. NEUTRINO SCATTERING 

In this section we review coherent neutrino scattering 
which is expected to dominate the neutrino opacity in 
regions where clusters, such as either nuclei or pasta, are 
present. Although the formalism has been presented al- 
ready in Ref. @, some details are repeated here (almost 
verbatim) for the sake of completeness and consistency. 

In the absence of corrections of order E v jM (with E v 
the neutrino energy and M the nucleon mass) and ne- 
glecting contributions from weak magnetism, the cross 
section for neutrino-nucleon elastic scattering in free 
space is given by the following expression |l5j : 



da 



^p[ C 2(3-cos0) + C 2(l + cos#)] , (7) 



where Gf is the Fermi coupling constant and 9 the scat- 
tering angle. 

Having neglected the contribution from weak mag- 
netism, the weak neutral current of a nucleon con- 
tains only axial-vector (757^) and vector 7^ contribu- 
tions. That is, 



The axial coupling constant is, 

Ca=±^ (9a = 1-26) 



(8) 



(9) 



Note that in the above equation the +(— ) sign is for 
neutrino-proton(neutrino-neutron) scattering. The weak 
charge of the proton c v is small, as it is strongly 
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suppressed by the weak-mixing (or Weinberg) angle 
sin 2 #vv = 0.231. It is given by 

c v = - - 2 sin 2 9 W = 0.038 « . (10) 

In contrast, the weak charge of a neutron is both large 
and insensitive to the weak-mixing angle: c„ = — 1/2. 

If nucleons cluster tightly, either into nuclei or into 
pasta, then the scattering of neutrinos from the various 
nucleons in the cluster may be coherent. As a result, the 
cross section will be significantly enhanced as it would 
scale with the square of the number of nucleons |16| . In 
reality, only the contribution from the vector current is 
expected to be coherent. This is due to the strong spin 
and isospin dependence of the axial current, which is ex- 
pected reduce its coherence. (Recall that in the nonrel- 
ativistic limit, the nucleon axial- vector current becomes 
l5l T z —> —GTz)- Since in nuclei and presumably also in 
the pasta most nucleons pair off into spin singlet states, 
their axial-vector coupling to neutrinos will be strongly 
reduced. Hence, in this work we focus exclusively on 
coherence effects from the vector current. Coherence is 
important in neutrino scattering from the pasta because 
the neutrino wavelength is comparable to the interparti- 
cle spacing and even to the intcrcluster separation. One 
must then calculate the relative phase for neutrino scat- 
tering from different nucleons and then add their con- 
tribution coherently. This procedure is embodied in the 
static structure factor S(q). 

The static structure factor per neutron is defined as 
follows: 

5 (q) = ^£|<*n|flq)l*o>| 2 , ( n ) 

where V&o and $ n are ground and excited nuclear states, 
respectively and the weak vector charge density is given 
by 

N 

P(q) = ^cxp(iq- Ti) . (12) 

i=l 

As the small weak charge of the proton [Eq. (|10|) ] will 
be neglected henceforth, the sum in Eq. Ill2[l runs only 
over the N neutrons in the system. The cross section per 
neutron for neutrino scattering from the whole system is 
now given by 

The above expression is the single neutrino-neutron cross 
section per neutron obtained from Eq. Q (with c a = 0) 
multiplied by S(q). This indicates that S(q) contains the 
effects from coherence. Finally, note that the momentum 
transfer is related to the scattering angle through the 
following equation: 

q 2 = 2El(l -cos 6) . (14) 



The static structure factor has important limits. For a 
detailed justification of these limits the reader is referred 
to Ref. 8]. In the limit that the momentum transfer to 
the system goes to zero (q — > 0) the weak charge density 
[Eq. I|12l) ] becomes the neutron number operator p(q = 
0) = N . In this limit the static structure factor reduces 
to, 

% = 0) = l((iV 2 )-(iV} 2 ). (15) 

Thus, the q— >0 limit of the static structure factor is re- 
lated to the fluctuations in the number of particles, or 
equivalently, to the density fluctuations. These fluctua- 
tions are themselves related to the compressibility and 
diverge at the critical point 01- As the density fluc- 
tuations diverge near the phase transition, the neutrino 
opacity may increase significantly. This could have a dra- 
matic effect on present models of stellar collapse. So far 
one dimensional simulations with the most sophisticated 
treatment of neutrino transport have not exploded [Tsj . 

In the opposite q — ► oo limit, the neutrino wavelength 
is much shorter than the interparticle separation and the 
neutrino resolves one nucleon at a time. This limit corre- 
sponds to quasielastic scattering where the cross section 
per nucleon in the medium is the same as in free space. 
Thus, the coherence disappears and 

S(q -> oo) = 1 . (16) 

In Monte Carlo as well as in molecular dynamics sim- 
ulations it is convenient to compute the static struc- 
ture factor from the neutron-neutron correlation func- 
tion g(r). Indeed, the static structure factor is obtained 
from the Fourier transform of the two-neutron correlation 
function. That is, 

5(q) = 1 + Pn J d?r(g{v) - l) exp(iq • r) . (17) 

The convenience of the two-neutron correlation function 
stems from the fact that it measures spatial correlations 
that may be easily measured during the simulations. It 
is defined as follows: 

1 N 

ff( r ) = T^E<*°W r - r «)l*°>' (18) 

where p n = N/V is the average neutron density. Opera- 
tionally, the correlation function is measured by pausing 
the simulation to compute the number of neutron pairs 
separated by a distance |r|. Note that the two-neutron 
correlation function is normalized to one at large dis- 
tances g(r — ► oo) = 1; this corresponds to the average 
density of the medium. 

IV. RESULTS 

In this section results are presented for a variety of 
neutron-rich matter observables over a wide range of den- 
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sities. Our goal is to understand the evolution of the sys- 
tem with density. From the low-density phase of isolated 
nuclei, through the complex pasta phase, to uniform mat- 
ter at high densities. All the results in this section have 
been obtained with an electron fraction and a tempera- 
ture fixed at Y e = 0.2 and T = 1 MeV, respectively. In 
core-collapse supernova the electron fraction starts near 
Y e — 0.5 and drops as electron capture proceeds. In a neu- 
tron star Y e is small — of the order of 0.1 — as determined 
by beta equilibrium and the nuclear symmetry energy (a 
stiff symmetry energy favors larger values for Y e |19|'). 
Thus, the value of Ye = 0.2 adopted here is representative 
of neutron-rich matter. 

There arc limitations in our simple semiclassical model 
at both low and high temperatures. At very low temper- 
atures the system will solidify while at high temperatures 
the model may not calculate accurately the free-energy 
difference between the liquid and the vapor. As a re- 
sult the pasta may melt at a somewhat too low of a 
temperature. Results are thus presented for only T = 1 
MeV where the model gives realistic results. Recall that 
our model reproduces both the saturation density and 
binding energy of nuclear matter, and the long-range 
Coulomb repulsion between clusters. Therefore, a good 
description of the clustering should be expected. 

The simulations start at the low baryon density of 
p = 0.01 fm -3 (which corresponds approximately to 
2 x 10 13 g/cm 3 ). This is about 1/15 of nuclear-matter 
saturation density. At this low density the most time 
consuming part of the simulation is preparing appropri- 
ate initial conditions. This is because the Coulomb bar- 
rier greatly hinders the motion of protons into and out 
of the clusters. The Coulomb barrier becomes an even 
greater challenge at lower densities. Thus, no effort has 
been made to simulate densities below 0.01 fm -3 . 

In Ref. results were presented for Monte Carlo sim- 
ulations with A = 4, 000 nucleons. Unfortunately, finite- 
size effects make it difficult to calculate S(q) accurately 
at small momentum transfers from these "small" simu- 
lations. Further, at a baryon density p = 0.01 fm -3 each 
side of the simulation volume has a length of approxi- 
mately L — 75 fm. This is inadequate as the wavelength 
of a 10 MeV neutrino is close to 120 fm. Therefore, in 
the present simulations the number of particles has been 
increased by a full of order of magnitude to A = 40, 000 
nucleons. This results in a simulation volume that has 
increased to a cube of L = 158.7 fm on a side. The 
molecular dynamics simulation starts with the nucleons 
uniformly distributed throughout the simulation volume 
and with a velocity profile corresponding toaT=l MeV 
Boltzmann distribution. The system is then evolved ac- 
cording to velocity- Verlet algorithm using a time step of 
the order of At — 2 fm/c. Velocity- Verlet is an energy 
conserving algorithm and with this time step the total 
energy of the system is conserved to one part in 10 5 . 
However, in order to preserve the temperature fixed at 
T=l MeV, the velocities of all the particles must be con- 
tinuously rescaled so that the kinetic energy per particle 



stays approximately fixed at ik-sT/2. For some excellent 
references to molecular dynamics simulations we refer the 
reader to E3, E3, E3 - 

In an attempt to speed up equilibration, the temper- 
ature of the system was occasionally raised during the 
evolution to 1.5 — 2 MeV. This could aid the system move 
away from local minimum, as is conventionally done with 
simulated annealing. Equilibration is checked by moni- 
toring the time dependence of the two-neutron correla- 
tion function g(r) and the static structure factor S(q). 
The peak in S(q) was observed to grow slowly with time 
as the cluster size increased. Changes became exceed- 
ingly slow after evolving the system for a long total time 
of 1,287,000 fm/c. Although long, further changes in 
S(q) over much larger time scales can not be ruled out. 
This suggests that our p = 0.01 fm -3 results may in- 
clude a systematic error due to the long — but finite — 
equilibration time. Fortunately, equilibrium seems to be 
reached much faster at higher densities so that the slow 
equilibration probably ceases to be a problem at these 
densities. 

Results for two observables — the potential energy per 
particle and the pressure — are displayed in Table [H] as a 
function of density. Also included in the table are the 
number of nucleons and the total evolution time for each 
density. Note that the pressure of the system is computed 
from the virial equation as follows [l2|: 



P 



_L/ v — 

3A\^ ni dr 



(19) 



In the case of non-interacting particles, Eq. Ijl9(l reduces 
to the well-known equation of state of a classical ideal gas. 
Thus, the second term in the above equation reflects the 
modifications to the ideal-gas law due to the interactions. 



TABLE II: molecular dynamics simulation results. Here p (in 
fm -3 ) is the baryon density, A is the baryon number, tf (in 
fm/c) is the total evolution time, V/A (in MeV) is the poten- 
tial energy per particle, P (in MeV/fm 3 ) is the pressure, and 
S(0) is the approximate value of the static structure factor at 
q — Q computed as in Eq. I26H . 



p 


A 


t f 


V/A 


P 


5(0) 


0.010 


40, 000 


1,287, 000 


-5.377(1) 


6.9 xlO" 3 


0.790 


0.025 


100, 000 


52,000 


-5.145(1) 


3.2 xlO" 2 


0.344 


0.050 


100, 000 


28, 000 


-4.463(1) 


0.13 


0.139 


0.075 


40, 000 


60, 000 


-3.686(1) 


0.33 


0.077 



The neutron-neutron correlation function g(r) at a 
baryon density of p = 0.01 fm -3 is shown in Fig.^ The 
two-neutron correlation function measures the probabil- 
ity of finding a pair of neutrons separated by a fixed dis- 
tance r. A large broad peak is observed in g(r) in the 2- 
10 fm region; the lack of neutrons with a relative distance 
of less than 2 fm is due to the hard core of the potential. 
The sharper sub-peaks contained in this structure reflect 



6 



BO 





1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 








[J — U.U1 1111 








Ye=0.2 








T=l MeV 












1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



10 20 30 40 50 60 70 80 
r(fm) 

FIG. 1: Neutron-neutron correlation function g(r) at a den- 
sity of p = 0.01 fm~ 3 , an electron fraction of Y e — 0.2, and 
a temperature of T = 1 MeV. This from a simulation with 
40,000 nucleons. 



neutron-neutron correlations (nearest neighbors, next-to- 
nearest neighbors, and so on) within a single cluster. The 
Coulomb repulsion among protons prevents the clusters 
from growing arbitrarily large and keeps them apart. The 
dip in g(r) at r ~ 10 fm is a result of the Coulomb re- 
pulsion between clusters. Finally, the small broad peaks 
near 25-30, 50, and 75 fm reflect correlations among the 
different clusters. 



Figure [21 displays an enlargement of the neutron- 
neutron correlation function for large values of r. Finite- 
size effects lead to an abrupt drop in g(r) at r = L/2~ 
80 fm (not shown). To ensure a reliable estimate of 
its Fourier transform — and correspondingly of the static 
structure factor S(q) — one must extrapolate g(r) to the 
region r > L/2. To do so, an analytic function of the 
following form is fitted to g(r): 

g a {r) = A e~ aor cos(k r + S ) + 1 . (20) 

The constants Aq, ao, fco, and So are obtained from a fit 
to the large-r behavior of the neutron-neutron correlation 
function. The result of this fit is indicated by the red 
dashed line in Fig. El 




FIG. 3: (Color online) The 0.03 fm" 3 proton density lsosur- 
face for one configuration of 40,000 nucleons at a density of 
0.01 fm -3 . The simulation volume is a cube 159 fm on a side. 
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g(r)=l forr>L/2 
g(r)=g () (r) for r>L/2 



0.99 
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50 60 70 80 90 100 

r(fm) 



FIG. 2: (Color online) The large-r behavior of the two- 
neutron correlation function displayed in Fig. Q Note the 
expanded y-scale. Also shown (red dashed line) is the ana- 
lytic fit to g(r) according to Eq. l20> . 



Figure |21 shows the 0.03 fm~ 3 isosurface of the proton 
density for one configuration of 40,000 nucleons at a den- 
sity of 0.01 fm -3 . All of the protons and most of the 
neutrons are clustered into nuclei. We discuss the size 
of these nuclei in Section El There is also a low density 
neutron gas between the clusters which is not shown. 

The static structure factor S(q) may now be calculated 
from the Fourier transform of g(r) [see Eq. (|17|) ]. That 



is, 



S(q) = 1 + Anp 7 , 



sin(fjrr) 



qr 



(g(r) - l)r 2 dr . (21) 



In Fig.^Jthe static structure factor obtained by using the 
above extrapolation (i.e., with go(r) for r>L/2) is dis- 
played with the red dashed line. In an earlier publication 
the static structure factor was calculated directly from 
the simulation results assuming g(r) = l for r>L/2 0. 
This result is also shown for comparison (black solid line) . 
The improvement in the low-momentum transfer behav- 
ior of S(q) is clearly evident. This is important as the 
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value of the static structure factor at zero-momentum 
transfer S(q= 0) monitors density fluctuations in the sys- 
tem and these may be indicative of a phase transition. 
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FIG. 4: (Color online) Static structure factor S(q) at a density 
of p = 0.01 fm~ 3 , an electron fraction of Y e = 0.2, and a 
temperature of T= 1 MeV. The black solid line assumes g(r) = 
1 for r > L/2 while the red dashed line includes an analytic 
extrapolation for g(r) for r>L/2. The green solid square is 
the value of 5(0) from Eq. iggj and Table El 



with the pressure expressed in units of MeV/fm 3 and the 
density in fm -3 . This yields, 



S(q = 0) 



1.358x10" 



(1.037 x lO- 3 + p a ) 



(26) 



These approximate values for S(q = 0) have been reported 
in Table ITT1 For comparison, they have also been added 
(with a green solid square) to the various figures for which 
S(q) was directly computed from the Fourier transform 
of the two-neutron correlation function (see Figs. |H| 
and !12|l . Note that there is good agreement between the 
two approaches. 

At low momentum transfers the static structure factor 
is small because of ion screening. Coulomb correlations — 
which both hinder the growth of clusters and keep them 
well separated — screen the weak charge of the clusters 
thereby reducing S(q). Further, a large peak is seen in 
S(q) for g ~ 0.25 fin -1 . This corresponds to the coher- 
ent scattering from all of the neutrons in a cluster. As 
we will show in the next section, this peak reproduces 
coherent neutrino-nucleus elastic scattering at low den- 
sities. Finally, the static structure factor decreases for 
q > 0.3 fm" 1 . This is due to the form factor of a clus- 
ter. As the momentum-transfer increases, the neutrino 
can no longer scatter coherently form all the neutrons 
because of the size of the cluster is larger than the neu- 
trino wavelength. All these features will be discussed in 
greater detail in the next section. 



In the limit of zero-momentum transfer the static 
structure factor S(q = 0) is directly related to the isother- 
mal compressibility. That is [l7j . 



S(q = 0) = /ofceT K-t = k?,T 



'dP^ 
dp / T 

where the isothermal compressibility is given by 

'dP\ (dp_ s 



JCrp 1 — 



-V 



(22) 



(23) 



For a classical ideal gas (i.e., P = pk-^T) the isothermal 
compressibility reduces to /C^ 1 = pk^T and S(q = Q) = 1. 
As expected, in the absence of interactions there are no 
spatial correlations among the particles. Assuming now 
that as q — > the fluctuations in the neutron density 
are proportional to the corresponding fluctuations in the 
baryon density, one obtains for the static structure factor 
per neutron 



S(q = 0) 



N. 
A' 



-knT 



dP 
~dp~ 



(24) 



The derivative of the pressure with respect to the baryon 
density has not been directly calculated in the simula- 
tions. However, the pressure has been computed at vari- 
ous densities and has been tabulated in Table ITT1 These 
values can be approximated by a simple fit of the form: 

P{p) = (0.611)/?+ (228. 131)p Q+1 (a = 1.583), (25) 
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FIG. 5: Neutron-neutron correlation function g(r) at a den- 
sity of p — 0.025 fm~ 3 , an electron fraction of Y e = 0.2, and 
a temperature of T — 1 MeV. This from a simulation with 
100,000 nucleons. 



Next, simulation results are presented at the higher 
density of p = 0.025 fm~ 3 . At this density it becomes 
much easier to equilibrate the system as protons have 
shorter distances to move over the Coulomb barriers. 
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FIG. 6: (Color online) The large-r behavior of the two- 
neutron correlation function displayed in Fig. |S] Also shown 
(red dashed line) is the analytic fit to g(r). 



In order to minimize finite-size effects more nucleons — 
a total of A = 100, 000 — are used for these simulations. 
The two-neutron correlation function is shown in Fig. |5j 
with its behavior at large distances amplified in Fig. 
Note that the calculation for g(r), which proceeds by 
histograming distances between the N(N—l)/2 pairs of 
neutrons, is now considerably more time consuming. The 
resulting static structure factor is shown in Fig. |SJ The 
peak in S(q) has now moved to larger q because of the 
shorter distance between clusters. 




FIG. 7: (Color online) The 0.03 fm" 3 proton density lsosur- 
face for one configuration of 100,000 nucleons at a density of 
0.025 fm~ 3 . The simulation volume is a cube 159 fm on a 
side. 



Figure shows the 0.03 fm~ 3 isosurface of the proton 
density for one configuration of 100,000 nucleons at a 
density of 0.025 fm~ 3 . All of the protons and most of the 
neutrons are clustered into nuclei. The size of these nuclei 
is now larger than at a density of 0.01 fm -3 as discussed 
in Section There is also a low density neutron gas 
between the clusters which is not shown. 
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FIG. 8: (Color online) Static structure factor S(q) at a density 
of p = 0.025 fm -3 , an electron fraction of Y e — 0.2, and a 
temperature of T — 1 MeV. The black solid line assumes g(r) = 
1 for r>L/2 while the red dashed line includes an analytic 
extrapolation for g(r) for r>L/2. The green solid square is 
the value of S(0) from Eq. and Tabic |Tj] 



Simulations have also been performed at a density of 
p = 0.05 fm~ 3 using A = 100, 000 nucleons. The two- 
neutron correlation function, together with its amplifi- 
cation at large values of r, are shown Figs. and 1101 
respectively. The corresponding static structure factor 
is displayed Fig. ^| Note the significant improvement 
in the behavior of S(q) at low momentum transfers as 
the sharp cutoff in g(r) is removed in favor of a smooth 
extrapolation [see Eq. l(2Tj|) ]. As the density increases, 
and thus the separation between clusters decreases, the 
peak in S(q) continues to move to higher q. However, the 
peak value of S(q) has now been reduced because of the 
increase in ion screening with density. 

Figure ITU shows the 0.03 fm -3 isosurface of the pro- 
ton density for one configuration of 100,000 nucleons at a 
density of 0.05 fm -3 . The clusters are now seen to have 
very elongated shapes. The low density neutron gas be- 
tween these clusters is not shown. 

We conclude this section by presenting results for the 
two-neutron correlation function at a density of p = 
0.075 fm~ 3 using a total of A = 40, 000 nucleons in 
Fig. ^3 Note that this is the largest density consid- 
ered in this work. At this density the clusters have 
been "melted" and the system has evolved into a uniform 
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FIG. 9: Neutron-neutron correlation function g(r) at a den- 
sity of p = 0.05 fm -3 , an electron fraction of Y e =0.2, and a 
temperature of T = 1 MeV. This is from a simulation with 
100,000 nucleons. 




r(fm) 

FIG. 10: (Color online) The large-r behavior of the two- 
neutron correlation function displayed in Fig. |5] Also shown 
(red dashed line) is the analytic fit to g(r). 



phase. In Fig.^Jresults for the static structure factor at 
this density are compared with the corresponding results 
at the lower densities. There is no longer evidence for 
a large peak in S(q) in the uniform system as a conse- 
quence of the complete loss in coherence. Finally, Fig. ^| 
shows S(q) at large momentum transfers. One observes 
that the static structure factor decreases with increasing 
density in the intermediate q-region before approaching 
the value of one (as it must) at high q. 




FIG. 11: (Color online) The 0.03 fm~ 3 proton density isosur- 
face for one configuration of 100,000 nucleons at a density of 
0.05 fm~ 3 . The simulation volume is a cube 126 fm on a side. 
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FIG. 12: (Color online) Static structure factor S(q) at a den- 
sity of p — 0.05 fm -3 , an electron fraction of Y e — 0.2, and 
a temperature of T = 1 MeV. The black solid line assumes 
g(r) — 1 for r > L/2 while the red dashed line includes an 
analytic extrapolation for g(r) for r>L/2. The green solid 
square is the value of S(0) from Eq. ^ and Table El 



In the next section we will compare these "complete" 
results with conventional approaches that model the sys- 
tem as a collection of strongly-correlated, neutron-rich 
nuclei plus a neutron gas. In this approach the peak 
observed in the static structure factor is attributed to 
neutrino-nucleus elastic scattering. Note, however, that 
the complete simulation results obtained in this section 
should remain valid even when these nuclear models 
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FIG. 13: Neutron-neutron correlation function g(r) at a den- 
sity of p — 0.075 fm~ 3 , an electron fraction of Y e = 0.2, and 
a temperature of T = 1 MeV. This from a simulation with 
40,000 nucleons. 




FIG. 14: (Color online) Static structure factor S(q) for a va- 
riety of densities at an electron fraction of Y e = 0.2 and a tem- 
perature of T=l MeV. The black solid line is for a density of 
p = 0.01 fm -3 , while the red dashed line is for p = 0.025 fm -3 , 
the green dotted line for p — 0.05 fm -3 , and the blue dot- 
dashed line for p = 0.075 fm" 3 . 



break down. 
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FIG. 15: (Color online) The large-g behavior of the static 
structure function displayed in Fig. 1141 



V. CLUSTER MODEL 

In this section a cluster model is developed with the 
goal of comparing our simulation results from the pre- 
vious section with commonly used approaches. In the 
complete model employed earlier, trajectories for all the 
nucleons were calculated from molecular dynamics simu- 
lations and these were used to compute directly the two- 
neutron correlation function and its corresponding static 
structure factor. One of the main virtues of such an 
approach is that there is no need to decide if a given nu- 
cleon is part of a cluster or part of the background vapor. 
Nevertheless, in this section a clustering algorithm is con- 
structed with the aim of assigning nucleons to clusters. 
In this way one can compare the inferred composition ex- 
tracted from our simulations with many nuclear statisti- 
cal equilibrium (NSE) models that describe the system as 
a collection of nuclei and free nucleons. In this way one 
can then compare the static structure factor extracted 
from the complete simulations with that calculated in 
these NSE models. 



A. Clustering Algorithm 

The clustering algorithm implemented in this section 
assigns a nucleon to a cluster if it is within a distance Rc 
of at least one other nucleon in the cluster. In practice, 
one stars with a given nucleon and searches for all of its 
"neighbors" , namely, all other nucleons contained within 
a sphere of radius Rc- Next, one repeats the same pro- 
cedure for all of its neighbors until no new neighbors are 
found. This procedure divides a fixed configuration of 
nucleons into a collection of various mass clusters (i.e., 
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TABLE III: Distribution of cluster as a function of the cutoff 
radius Rc at a density p = 0.01 fm~ 3 for a system of 40,000 
nucleons. The number of free neutrons is denoted by ^(^4 = 
1), the average mass of all A>2 clusters by (A), and the size 
of the largest cluster by v4 m ax. 



Rc (fm) 


N(A = 1) 


(A) 




2.0 


20,585 


46.20 


118 


2.5 


14,233 


102.29 


155 


3.0 


11,062 


98.75 


160 


3.5 


7,856 


96.17 


254 


4.0 


5,131 


100.97 


261 


4.5 


2,888 


149.84 


513 


5.0 


1,427 


14,069 


23,189 



"nuclei"). 

To illustrate this procedure the final nucleoli configura- 
tion of the complete simulation of the previous section at 
a density of p = 0.01 fm -3 is selected after the system has 
evolved for a total time of tf = 1, 287, 000 fm/c. Having 
selected a cutoff radius of Rc = 3 fm, one finds that the 
40,000 nucleons in the system are divided in the follow- 
ing way: a) 11,062 free neutrons, b) no free protons, c) a 
few light nuclei with A<8, and d) a collections of heavy 
nuclei with a mass distribution of 50 < A < 160. The 
mass- weighted average of all clusters with A > 2 is equal 
to (A) = 99. This distribution of clusters is displayed in 
Fig. [inland listed in Table ITTT1 
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TABLE IV: Distribution of cluster as a function of the cutoff 
radius Rc at a density p = 0.025 fm -3 for a system of 100,000 
nucleons. The number of free neutrons is denoted by ^(^4 = 
1), the average mass of all A>2 clusters by {A), and the size 
of the largest cluster by A max . 



Rc (fm) 


N(A = 1) 


(A) 




2.0 


49,761 


58.95 


162 


2.5 


28,234 


167.40 


423 


3.0 


14,549 


198.91 


614 


3.5 


5,187 


59,411 


75,003 


TABLE V: Distribution of cluster as a function of the cutoff 
radius Rc at a density p = 0.05 fm -3 for a system of 100,000 
nucleons. The number of free neutrons is denoted by ^(^4 = 
1), the average mass of all yl>2 clusters by (A), and the size 
of the largest cluster by 4 max . 


Rc (fm) 


N(A = 1) 


(A) 


^4 max 


2.0 


47,239 


56.86 


308 


2.5 


16,219 


72,952 


78,178 



this range should give similar results. However, if Rc is 
chosen too large, for example Rc = 5 fm, then most of 
the nucleons become part of one single giant cluster. 
Similar results for a density of p = 0.025 fm -3 are pre- 
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sented in Table lIVl This distribution is extracted from 
the final configuration of 100,000 nucleons obtained after 
a total evolution time of tf = 52, 000 fm/c. Using a cutoff 
radius of Rc = 3 fm, the 100,000 nucleons in the system 
are now divided into 14,549 free neutrons, no free pro- 
tons, a few light nuclei, and a broad collection of heavy 
nuclei with A from about 80 to 614 nucleons. Such a 
distribution is shown in Fig. 1171 The average mass has 
now grown to (A) = 199. The mass of the heavy nuclei 
is seen to increase with density as shown in Fig. El Fi- 
nally, results at p = 0.05 fm -3 are presented in Table [V] 
Now the density is so high that it is difficult to design 
a sensibly scheme to divide the system into clusters, see 
Fig. ^] For example, even with a cutoff radius as small 
as R c = 2.5 fm, already 78,178 of the 100,000 nucleons 
become part of a single giant cluster. 



Cluster Form Factors 



FIG. 16: Number of clusters of atomic mass A for one config- 
uration of 40,000 nucleons at a density of p = 0.01 fm~ 3 . Note 
that this is a linear- log plot. 



In Table ITTT1 results are also displayed for values of the 
cutoff radius in the 2-5 fm range. Note that the average 
cluster mass (A) appears remarkably constant for 2.5 < 
Rc < 4 fm. This suggest that any value of Rc within 



To describe coherent neutrino scattering from a single 
cluster, that is, neutrino-nucleus elastic scattering, one 
must calculate the elastic form factor for the cluster. This 
is given by 



JV 



Hi) 



1 Sin (9 r ") 



N 



(27) 



Here the sum runs over the N neutrons in the cluster and 
r n is distance from the n t h neutron to the center of mass 
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FIG. 17: (Color online) Number of clusters of atomic mass A 
for one configuration of 100,000 nucleons at a density of p = 
0.025 fm~ 3 (black hatched line). Also shown for comparison 
is the number of clusters (from Fig, lib! at p = 0.01 fm -3 (red 
solid line). Note that both scales are now linear. 



the exotic, neutron-rich nucleus 98 Ni calculated in a rela- 
tivistic mean-field a ppr oximation using the very success- 
ful NL3 interaction [2£|. While the NL3 form factor has 
a slightly smaller RMS radius, the overall agreement be- 
tween all three models is fairly good. Note that 98 Ni 
(rather than 100 Ni) was used in this calculation as it 
contains closed protons and neutrons subshells. For this 
exotic nucleus the l/i 11//2 neutron orbit — responsible for 
magic number 82 — is not even bound. 




of the TV-neutron system. The form factor represents the 
Fourier transform of the point neutron density and here, 
for simplicity, has been averaged over the direction of the 
momentum transfer. Note that the elastic form factor is 
normalized so that F(q = 0) = 1. In Fig. ^|the elastic 
form factors of all clusters with A > 10 are displayed at 
a density of p = 0.01 fm -3 (a cutoff radius of Rc = 3 fm 
was selected). The large spread in the form factors re- 
flects the many different sizes of the individual clusters 
(see Fig. I16|l . Indeed, the root-mean-square (RMS) ra- 
dius of a cluster appears to scale approximately as A 1 / 3 . 
Therefore, in Fig. ^j] all of these form factors are plot- 
ted but against a scaled momentum transfer q A 1 / 3 . Now 
all the (scaled) form factors fall in a fairly narrow band 
suggesting that, while these neutron-rich clusters have 
different radii, they all share a similar shape. 

In Fig. [20| we display the form factor for a single 
neutron-rich cluster with A = 100 and Z = 28 (" 100 Ni") 
extracted from the simulation with a density of p=0.01 
fm~ 3 . Also shown in the figure (with a green dotted line) 
is the form factor Fo(q) of a uniform neutron distribution 
with a sharp surface radius R n chosen to reproduce the 
RMS radius of the neutron-rich cluster < r\ > 1//2 . It is 
given by 

_ . , „sin(a;) — a;cos(a;) . . . 

F (q) = 3 y > ^ K -L (x = qR n ) , (28) 

with R n given by 

Rn = \j\irl) 1 " ■ (29) 
Finally, Fig. [201 also shows the neutron form factor of 



FIG. 18: Cluster form factor F(q) as a function of the mo- 
mentum transfer q for all clusters with A > 10 using a single 
configuration of 40,000 nucleons at a density of p — 0.01 fm -3 . 
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FIG. 19: Cluster form factor F(q) as a function of the scaled 
momentum transfer qA 1 ^ 3 for all clusters with A > 10 using a 
single configuration of 40,000 nucleons at a density of p = 0.01 
fm" 3 . 



The clusters generated in the simulations are neutron- 
rich nuclei with well developed neutron skins. Nuclei with 
neutron skins are characterized by neutron radii that are 
larger than those for the protons. Using the distribution 
of nuclei obtained with a density of p — 0.01 fm -3 , the 
following values are obtained for average matter, proton, 
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FIG. 20: (Color online) Cluster form factor F(q) as a function 
of the momentum transfer q. The red dashed line represents 
the angle-averaged form factor for one cluster with A — 100 
and Z = 28 from a simulation at a density of p = 0.01 fm~ 3 . 
The dotted (green) line is the form factor of a uniform density 
sphere with the same root-mean-square radius [see Eqs. 12811 
and 129H 1. Finally, the solid line is the form factor of the very 
neutron-rich nucleus 98 Ni calculated in a relativistic mean- 
field approximation with the NL3 interaction |20| . 



and neutron RMS radii respectively: 

<r 2 ) 1/2 = 1.06 A 1 / 3 fm . 
( r 2)V2 = 0.91 A 1 ' 3 fm : 
(r 2 ) 1 / 2 = 1.11 A 1 / 3 fm 



(30a) 
(30b) 
(30c) 



Note that the sharp-surface radius of a uniform distri- 
bution with the same RMS radius is simply given by 
(5/3) 1/2 times these values [see Eq. <|29"|) ]. 



C. Single Heavy Nucleus Approximation 

A number of approaches to dense matter, such as those 
using the equation of state by Lattimer and Swesty [2lj| . 
model the system as a collection of free neutrons plus 
a single representative heavy nucleus. Occasionally, free 
protons and alpha particles are also added to the sys- 
tem. To mimic this approach, a model is constructed 
based on our earlier cluster results reported in Tables ITTT1 
and IIVI for a cutoff radius of Rc — 3 fm. For example, 
at a density of ,0 = 0.01 fm -3 the system contains a mass 
fraction X n — 0.28 of free neutrons and a mass fraction 
of Xh = 1 — X n = 0.72 for the single representative heavy 
nucleus. According to the average mass reported in Ta- 
ble lllll a mass of A = 100 is assigned to this representative 
heavy nucleus. Conservation of charge constrains this nu- 
cleus to have Z 28. Note that due to the presence of free 



TABLE VI: Composition of the system in the single-heavy- 
nucleus approximation. The mass fraction of free neutrons 
is denoted by is X n and that of heavy nuclei by Xh- The 
mass and charge of the nuclei are given by A and Z, respec- 
tively. Finally, the radii of the equivalent uniform proton and 
neutron distributions are denoted by R p and R n , respectively. 



P (fm" 3 ) 


X n 


X h 


A 


Z R p (fm) 


R n (fm) 


0.010 


0.28 


0.72 


100 


28 5.45 


6.68 


0.025 


0.14 


0.86 


199 


47 6.84 


8.40 



neutrons (but not free protons) the charge-to-mass ratio 
of the heavy nucleus Z / A = 0.2& slightly exceeds the elec- 
tron fraction Y" e = 0.2 of the whole system. The assumed 
composition of the system at densities of p = 0.01 fm~ 3 
and p = 0.025 fm~ 3 is given in Table WH 

The heavy nuclei are assumed to interact exclusively 
via a screened Coulomb interaction. Each nucleus is as- 
sumed to have a uniform charge distribution p c h that 
extends out to a radius R p chosen to reproduce the pro- 
ton RMS radius (r 2 ) 1 / 2 given in Eq. (JSOJl. The Coulomb 
interaction between two such nuclei whose centers are 
separated by a distance R is given by 



V C (R) = e 2 / d 3 rp ch (r) / <fr' Pch {r') 



z—Rtot I A 
-Rtot 



(31) 



where i?t. t = |R + r — r'| and A is the screening length 
fixed (as in Sec. [H| at a constant value of A = 10 fm. In 
the limit that the distance between nuclei is much larger 
than the nuclear RMS radius (i.e., R p <C R) the above 
integral reduces to 



V C (R) c e 2 p 2 ch 



Z 2 e 2 



-R/X 

2 ,2 ' / / d 3 r e ~rcos8/\ 



R 



R 

,- R /*f{R p /\) 



(32) 



where the dimcnsionless function f(x) has been defined 
as follows: 



/(*) = 



l xcosh(a;) — sinh(a;) 



(x = R p /X) . (33) 



Note that the function / is independent of R. Indeed, it 
only depends on the dimcnsionless ratio R p /X, namely, on 
the interplay between the nuclear size and the screening 
length. In the absence of screening, / = 1 (independent 
of nuclear size) in accordance with Gauss' law. However, 
with screening / becomes greater than one. The finite 
nuclear size places some of the charges closer together 
than i?; this increases the repulsion. Of course, the fi- 
nite size also places some of the charges farther apart, 
thereby decreasing the repulsion. When these two effects 
are weighted by the screening factor e~ r / A , the repulsion 
more than compensates for the "attraction" leading ul- 
timately to / > 1. In the particular case of R p — 5.45 fm 
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and A=10 fm, one obtains / (0.545) = 1.061 (about a 6% 
increase). 

The single-heavy-nucleus models consists of a gas of 
nonintcracting neutrons plus ions interacting via the 
screened Coulomb interaction given in Eqs. (S2J}. Molec- 
ular dynamics simulations in the ion coordinates are per- 
formed to compute its static structure factor Si on (q). The 
simulations used 5,000 to 10,000 ions and a time step of 
25 to 75 fm/c. The ion simulation can afford a larger time 
step than the corresponding nucleon simulation because 
the heavier ions move slower. Further, the ion simula- 
tions require fewer particles to simulate the same physi- 
cal volume because each ion "contains" several nucleons. 
The static structure factor for the ions computed in this 
way (for densities of p = 0.01 fm -3 and p — 0.025 fm -3 ) 
is shown in Fig. 
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FIG. 21: (Color online) Ion static structure factor Si on (q) 
as a function of the momentum transfer q. The black solid 
line is from of a simulation with 10,000 ions corresponding 
to a density of p — 0.025 fm -3 . The red dashed line is from 
a simulation with 5,000 ions corresponding to a density of 
p = 0.01. See text for details. 

Neutrino scattering from this system is described, 
in this single-heavy-nucleus model, by neutrino-nucleus 
elastic scattering within a framework that incorporates 
effects from both, the nuclear form factor and ion screen- 
ing from the correlated nuclei. The cross section for elas- 
tic neutrino scattering from a single nucleus is propor- 
tional to the square of the weak charge of the nucleus 
times a suitable form factor to account for it finite size. 
For the weak charge of the nucleus we simply use its neu- 
tron number N as we continue to ignore the small weak 
charge of the proton, i.e., <2 WC ak = — N+Z(l— 4 sin 2 9w) — > 
— N. Thus, the weak nuclear form factor reduces to that 
of the neutron distribution. Further, to incorporate ef- 
fects that result from correlations among the ions, such 
as ion screening, the cross section is multiplied by the 



ion static structure factor Si on (q). Finally, one multi- 
plies these terms by the fraction Xh of heavy nuclei and 
divides by N to obtain a static structure factor per neu- 
tron S mo< i e \{g) consistent with the normalization of the 
earlier sections. That is, 

Smodcifa) = X h NF(q) 2 S ion {q) . (34) 

Note that in addition to the coherent nuclear contribu- 
tion there is a small incoherent contribution from the 
neutron gas that has been neglected. As defined above, 
this static structure factor can now be directly com- 
pared to the one obtained in the full nucleon simula- 
tions. This prescription for S mo dei(q) corresponds to what 
is presently used in most supernova simulations. These 
simulations often take Xh and N from the Lattimer- 
Swesty equation of state |2l| and Si on (q) as computed 
in Ref. H|. 

In Fig.|221thc model static structure factor S mo dei(<z) is 
compared to the one from the full nucleon calculation (see 
Sec. liVf) at a density of p = 0.01 fm -3 . The uniform form 
factor of Eq. (|28|l is used with the sharp surface radius 
R n listed in Table IVll For low to moderate momentum 
transfers the agreement between the two approaches is 
excellent. This indicates that — at this density and (low) 
momentum transfers — the system is well described by a 
collection of nuclei of a single average mass. We expect 
that this good agreement will also hold at lower densi- 
ties. However, there is a modest disagreement between 
Smodci^) and the complete S(q) for g>0.25 fm -1 . This 
provides the first indication of limitations within the sin- 
gle heavy nucleus approximation. The discrepancy could 
arise because the broad distribution of cluster sizes dis- 
played in Fig. 1161 is approximated by a single average 
cluster with a mass of A = 100. Or it could be due to 
a breakdown in the factorization scheme. That is, the 
cross section may no longer factor into a product of a 
correlation function between ions (Si on ) times the weak 
response of a single ion (N 2 F(q) 2 ). 

A similar comparison is done in Fig. l23l but now at the 
higher density of p = 0.025 fm~ 3 . Now the disagreement 
between 5 mo dci(<7) and S(q) is more severe. This indi- 
cates that errors in the single nucleus approximation will 
grow rapidly with density. Moreover, the single nucleus 
approximation overpredicts the neutrino opacity relative 
to the complete calculations. 

VI. CONCLUSIONS 

Nonuniform neutron-rich matter was studied via semi- 
classical simulations with an interaction that reproduces 
the saturation density and binding energy of nuclear mat- 
ter and incorporates the long-range Coulomb repulsion 
between protons. Simulations with a large number of 
nucleons (40,000 to 100,000) enable the reliable deter- 
mination of the two-neutron correlation function and its 
Fourier transform — the static structure factor — even for 
low momentum transfers. The static structure factor 
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FIG. 22: (Color online) Neutron static structure factor S(q) 
as a function of the momentum transfer q at a density of 
p = 0.01 fm for the full calculation (black solid line). Also 
shown (red dashed line) is the prediction from the ion static 
structure factor in Fig. including the square of the cluster 
form factor, as explained in the text. 



FIG. 23: (Color online) Neutron static structure factor S(q) 
as a function of the momentum transfer q at a density of 
p = 0.025 fm~ 3 (black solid line). Also shown (red dashed 
line) is the prediction from the ion static structure factor in 
Fig. 1211 including the square of the cluster form factor, as 
explained in the text. 



S(q) describes coherent neutrino scattering that is ex- 
pected to dominate the neutrino opacity. At low mo- 
mentum transfer q the static structure factor is small 
because of ion screening; correlations between different 
clusters screen the weak charge. At intermediate mo- 
mentum transfers a large peak is developed in S(q) cor- 
responding to coherent scattering from all of the neutrons 
in a cluster. This peak moves to higher q and decreases 
in amplitude as the density of the system increases. 

In principle the neutrino opacity could be greatly in- 
creased by large density fluctuations. A simple first- 
order phase transition has a two-phase coexistence region 
where the pressure is independent of density. Large den- 
sity fluctuations in this region imply a very large value 
of the static structure factor at very small momentum 
transfers. Indeed, S(q = 0) is directly proportional to the 
density fluctuations in the system. Moreover, density 
fluctuations are also proportional to the isothermal com- 
pressibility. For consistency, the static structure factor 
S(q = 0) was computed in these two equivalent yet in- 
dependent ways, namely, as the Fourier transform of the 
two-neutron correlation function and as the derivative of 
the pressure with respect to the baryon density. While we 
find good agreement between these two schemes, no evi- 
dence is found in favor of a large enhancement in S(q = 0). 
We conclude that the system does not undergo a simple, 
single component first-order phase transition, so no large 
increase in the neutrino opacity was found. 

To compare our simulation results to more conven- 
tional approaches of wide use in supernova calculations a 
cluster model was introduced. A minimal spanning tree 



clustering algorithm was used to determine the composi- 
tion of the various clusters ( "nuclei" ) in the simulations 
(see for example Ref . [23| ) . To make contact with some of 
these conventional approaches, such as the single heavy 
nucleus approximation, the neutrino opacity was com- 
puted in a system modeled as a gas of free neutrons and 
a representative (i.e., average) single-species heavy nu- 
cleus. The neutrino opacity for such a system is domi- 
nated by elastic scattering from the heavy nucleus. The 
contribution from the single nucleus to the neutrino re- 
sponse is proportional to the square of its weak charge 
(assumed to be carried exclusively by the neutrons) and 
its elastic neutron form-factor, that accounts for its finite 
size. Further, Coulomb correlations among the differ- 
ent nuclei was incorporated through an ion static struc- 
ture factor to account for ion screening. Fairly good 
agreement is found between the single heavy nucleus ap- 
proximation and our complete simulations at low density 
and especially at small momentum transfers. However, 
starting at a density of approximately 10 13 g/cm 3 , we 
find a large disagreement between the two approaches 
that grows rapidly with increasing density. In partic- 
ular, our complete simulations yield neutrino opacities 
that are smaller than those in the single heavy nucleus 
approximation. Note that our full simulations yield accu- 
rate results even at the (high) densities where the single 
heavy nucleus approximation becomes invalid. We reiter- 
ate that the single heavy nucleus approximation is what 
is presently employed in most supernova simulations. 

Future work could include calculating the dynamical 
response of the system to study the transfer of energy 
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between neutrinos and matter. Note that a great virtue 
of molecular dynamics approaches combined with spe- 
cial purpose computers (such as what has been done 
here) is that dynamical information for systems with 
large number of particles may be readily obtained from 
time-dependent correlations. Particularly interesting is 
the low-energy part of the response which should be 
dominated by the so-called Pygmy resonances. These 
oscillations of the neutron skin of neutron-rich nuclei 
against the symmetric core should be efficiently excited 
by low-energy neutrinos. Another promising area for fu- 
ture research is the spin response of the system. A first 
step could involve including spin dependent forces in our 
model. The spin response is interesting because nucleons 
have large spin dependent couplings to neutrinos. 

APPENDIX: MDGRAPE 

To do the simulations we used a special purpose com- 
puter called the MDGRAPE-2. The MDGRAPE-2 is a 
single board which plugs into the PCI bus of a general 
purpose computer, and is designed for extremely fast cal- 
culation of forces and potentials in molecular dynamics 
simulations j^. It is the third generation of such hard- 
ware, which evolved from the work of J. Makino et al at 
the University of Tokyo on similar hardware called the 
GRAPE (for GRAvity PipE), for doing gravitational N- 
body problems |25| . In our case, we have two boards 
plugged into the PCI bus of one of the Power3+ nodes 
of Indiana University's IBM SP supercomputer. Each 
board is rated at a peak speed of 64 GigaFLOPS (float- 
ing point operations per second). The MDGRAPE-2 can 
compute any central potential of the form 

V{i,j) = b ij f{a ij {rl+e%)) (A.l) 

or the corresponding central force, which is of the same 
form, except multiplied by r. All three terms in (2) are 
of this form. In our case ey = 0, and 6y and ay are 
either scalars, or 2 x 2 matrices, corresponding to the 
two particle types proton and neutron. The boards are 
accessed via the M2 library, which the user links into his 
code. The library is very easy to use, and handles dis- 
tribution of work between the two MDGRAPE-2 boards 
without user intervention. The user defines f(x) by a 
function table of 1024 points, which the MDGRAPE-2 
interpolates via fourth degree polynomial interpolation. 
One can thus reproduce most physically realistic func- 
tions very accurately. Software is provided for construct- 
ing function tables, which are stored in files and loaded 
in during runtime. 

At each MD time step, one calls M2 subroutines to load 
in the function table and the scale factors or matrices ay 
and bij . One then calls a subroutine to load in the source 
particle coordinates, and subroutines to load integer ar- 
rays of particle types (0 for neutron, 1 for proton) for 
both source and target particles. Then one calls a force 
calculation routine, passing it the array of target particle 



coordinates. In our case the source and target particles 
are the same, but they do not have to be. One input 
parameter to the force calculation specifies that periodic 
boundary conditions should be used. The MDGRAPE-2 
has built in hardware for taking periodic b.c. into ac- 
count. The output is an array containing the total force 
on each target particle. A similar call can be made to 
compute the total potential energy of each particle. 

We must go through these steps three times, once for 
each term in (2). Still, the MDGRAPE-2 is much faster 
than serial Fortran code. In ordinary Fortran, this whole 
calculation would be done in a pair of nested DO loops, 
and thus would take of order 0{A 2 ) time. For our simula- 
tions with A = 40, 000 this would be prohibitive even for 
today's fast CPUs. But the two MDGRAPE-2 boards to- 
gether can do the force calculation about 90 times faster 
than a single Power3+ processor, so that a simulation of 
100,000 MD time steps that would take over two years 
using a serial program can be done in less than nine days. 
Benchmark tests show this speedup holds out to at least 
A = 160,000. Each MDGRAPE-2 board has enough 
memory to hold a half million particles, so we have not 
yet reached our maximum capability. 

We calculated the neutron-neutron correlation func- 
tion g(r) (see below) using ordinary Fortran code, as it 
was not clear how to perform this calculation with the 
MDGRAPE-2. Although g{r) is also an 0(A 2 ) calcula- 
tion, it is done infrequently, and does not severely impact 
performance. 
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